###################################################################
#Figure 2a - Birthplace and Representation Combination Predictions
###################################################################

library(foreign)
library(ggplot2)
library(lattice)
library(Hmisc)
library(readstata13)
library(latticeExtra)
library(gridExtra)
#loadfonts()
ssdata<-read.dta13("SameBirthSameRepRegion_ContOnly_Predictions.dta")
sddata<-read.dta13("SameBirthDiffRepRegion_ContOnly_Predictions.dta")
dsdata<-read.dta13("DiffBirthSameRepRegion_ContOnly_Predictions.dta")
dddata<-read.dta13("DiffBirthDiffRepRegion_ContOnly_Predictions.dta")
scaleD<-dddata[,1]
scaleS<-ssdata[,1]
LowBoundSS<-ssdata[,2]
UpBoundSS<-ssdata[,3]
LowBoundSD<-sddata[,2]
UpBoundSD<-sddata[,3]
LowBoundDS<-dsdata[,2]
UpBoundDS<-dsdata[,3]
LowBoundDD<-dddata[,2]
UpBoundDD<-dddata[,3]

xdiff<-1:2900
xsame<-1:1900

ydd<-append(LowBoundDD, rev(UpBoundDD))
yds<-append(LowBoundDS, rev(UpBoundDS))
ysd<-append(LowBoundSD, rev(UpBoundSD))
yss<-append(LowBoundSS, rev(UpBoundSS))
xmilesd<-append(scaleD, rev(scaleD))
xmiless<-append(scaleS, rev(scaleS))

xleg<-c(0, 3000, 3000, 0)
yleg<-c(.462, .462, .465, .465)

xlegbox<-c(100, 260, 260, 100)
ylegbox1<-c(0.4643, 0.4643, 0.4637, 0.4637)
ylegbox2<-c(0.4633, 0.4633, 0.4627, 0.4627)

masterD<-as.data.frame(cbind(xmilesd, ydd, yds))
masterS<-as.data.frame(cbind(xmiless, ysd, yss))
legendplot<-as.data.frame(cbind(yleg, xleg))
legendbox<-as.data.frame(cbind(xlegbox, ylegbox1, ylegbox2))

pdf(file="Predictions_LatticeContOnly_Figure2a.pdf", family="Bookman")
p<-ggplot(masterD, aes(x=xmilesd, y=ydd))+geom_polygon(fill="gray42")
p<-p+geom_polygon(data=masterD, aes(x=xmilesd, y=yds), fill="gray72")
p<-p+scale_y_continuous(name="Probability of Agreement", limits=c(.435, .465))+scale_x_continuous(name="Distance Between Birth Places (in miles)")+theme_bw()+theme(axis.text.x=element_text(angle=-90, vjust=.5), axis.title.x=element_text(size=9))+ggtitle("Different Birth Regions")+theme(plot.title=element_text(hjust=.5))
p<-p+geom_polygon(data=legendplot, aes(x=xleg, y=yleg), color="black", fill="white")
p<-p+geom_polygon(data=legendbox, aes(x=xlegbox, y=ylegbox1), fill="gray42")
p<-p+geom_polygon(data=legendbox, aes(x=xlegbox, y=ylegbox2), fill="gray72")
p<-p+annotate("text", x=1650, y=0.464, label="Different Regions of Representation", size=3)
p<-p+annotate("text", x=1500, y=0.463, label="Same Region of Representation", size=3)
p1<-ggplot(masterS, aes(x=xmiless, y=ysd))+geom_polygon(fill="gray42")
p1<-p1+geom_polygon(data=masterS, aes(x=xmiless, y=yss), fill="gray72")
p1<-p1+scale_y_continuous(name=" ", limits=c(.435, .465))+scale_x_continuous(name="Distance Between Birth Places (in miles)")+theme_bw()+theme(axis.text.x=element_text(angle=-90, vjust=.5), axis.title.x=element_text(size=9))+ggtitle("Same Birth Region")+theme(plot.title=element_text(hjust=.5))
pcomb<-grid.arrange(p, p1, nrow=1)
pcomb
dev.off()



############################################
#Figure 2b - Age Effects
############################################

library(foreign)
library(extrafont)
library(readstata13)
loadfonts()
data5<-read.dta13("Age5th_ContOnly_Predictions.dta")
dataM<-read.dta13("AgeMedian_ContOnly_Predictions.dta")
data95<-read.dta13("Age95th_ContOnly_Predictions.dta")
scale<-data5[,1]
LowBound5<-data5[,2]
UpBound5<-data5[,3]
LowBoundM<-dataM[,2]
UpBoundM<-dataM[,3]
LowBound95<-data95[,2]
UpBound95<-data95[,3]
pdf(file="Predictions_ContOnly_AgeEffs_Figure2b.pdf", family="Bookman")
plot(0, 1, type="n", xlim=c(0,2900), ylim=c(.444, .46), xlab="Distance Between Birth Places (in miles)", ylab="Probability of Agreement")
title(main=" ")
box() 

polygon(c(scale, rev(scale)), c(LowBound95, rev(UpBound95)), col="gray42", border=NA)

polygon(c(scale, rev(scale)), c(LowBound5, rev(UpBound5)), col="gray72", border=NA)

polygon(c(scale, rev(scale)), c(LowBoundM, rev(UpBoundM)), density=20, angle=45, col="black")

text(2310, .46, "Difference in Age", font=3)

legend("topright", cex=.98, legend=c("  ", "5th Percentile (1 year)", "Median (10 years)", "95th Percentile (28 years)"), fill=c("white", "gray72", "black", "gray42"), density=c(NA, NA, 30, NA), angle=c(NA, NA, 45, NA), border=c(NA, NA, "black", NA))

dev.off()

